Программа нахождения 94 цифр числа e.
Первый вариант вычисления цифр числа e указан в журнале Наука и Жизнь №4 1990 г.: вырезка в формате djvu. Там получилось только 27 цифр. Через год с небольшим опубликовали вариант, где цифр уже стало 92: Наука и жизнь №6 1991.
По сути, это уже максимум. При вычислениях (нахождение остатка от деления), чтобы не потерять точность при переполнении, приходится хранить только 6 цифр, а не 7. В результате получаем 15 регистров × 6 знаков = 90. Начало используется с переполнением, поэтому 91 цифра после запятой или всего 92.
Текущая версия – это улучшенный вариант того, что было опубликовано в 1991:
С учетом того, что в конце вычислений уже можно не боятся переполнений, в последний регистр дописаны ещё две цифры числа e. Так и получается 94 цифры.
Сначала поясню загадочное число 65, используемое в программе, но мало освещённое в исходной статье. Дело в том, что точность суммы сходящегося ряда, определяется последним членом. При вычислении e по формуле ряда ∑ 1⁄N!, точность определяется последним 1⁄N! Для определения разрядности (десятичной) любого натурального числа используется десятичный логарифм. Количество разрядов у числа определяется значением логарифма. Если учесть, что N! = 1×2×3×…×N, а логарифм произведения равен сумме логарифмов отдельных членов, то достаточно вычислить lg(1)+lg(2)+…+lg(N) для получения количества знаков у числа N! Например, если результат больше трех, то число содержит больше трех знаков (четыре).
Для вычислений суммы логарифмов составим простую программу:
| # | | 00 | 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 |
|---|---|---|---|---|---|---|---|---|---|---|
| 00 | | x→П0 | Cx | П→x0 | FLg | + | FL0 | 02 | С/П | ||
На входе число N, на выходе сумма логарифмов. Зададим 4 и получим 1.3802112, т.е. 2 знака, как и ожидалось: 1×2×3×4 = 24.
А если задать 65, то получим 90.916331, те самые 91 знак после запятой. Вот откуда это число 65 – это предел суммы ряда, чтобы получилась 91 точная цифра.
Здесь он раскрыт, потому что в исходной статье изложен кратко.
Используя ряд Тейлора, как сумму обратных факториалов, но разложенную по схеме Горнера: 1 + 1/2(1 + 1/3(1 + 1/4(…))), на каждом шаге вычисляем 1/N. Для дальнейшего понимания, как делим на N частями, используем пример. Пуcть хотим вычислить 1/73. Возьмём вместо единицы 108 и поделим группами по 2 знака (102=100).
| Делимое | Делитель | Целочисленный результат | Остаток |
|---|---|---|---|
| 100 | 73 | 01 | 27 |
| 2700 | 73 | 36 | 72 |
| 7200 | 73 | 98 | 46 |
| 4600 | 73 | 63 | 1 |
Таким образом, получим цифры результата (из 3-й колонки) 01369863 или с учётом заимствованных 108 это 1.3698630-02. Если вычислить 1/73 напрямую, то получим то же самое.
В нашем случае вместо 108 будем использовать 1091.
Делать будем целочисленно группами по 6 цифр, первая группа будет содержать 7.
Получится 107 и 14 раз по 106.
В каждой группе результат (целый, 6 знаков) сохраняется в регистре,
а остаток переносится
и умножается на 106.
Таким образом пробегаем по всем группам (регистрам). Итог деления - в регистрах, и остаток в стеке.
После деления на всё разряды остаток отбрасывается – не влияет уже.
Повторяем N (65) раз с уменьшением N.
Алгоритм целочисленного деления и нахождения остатка выполняется в подпрограмме. Накопленный остаток, как и N, хранится в стеке, так же как и в исходной программе. Основную часть исходной программы занимает многократный повтор команд: извлечь из очередного регистра текущие разряды, вызвать подпрограмму, сохранить. После пробега по всем разрядам уменьшить число N. При обнулении – закончить.
Чтобы сократить длину программы, многократный вызов подпрограммы с разными регистрами сделан во внутреннем цикле. А для счётчика этого цикла, как и возможность косвенной адресации номера регистра, использовано два последних разряда регистра Re, потому что при вычислениях значимы только первые шесть. За счёт этого программа значительно сократилась.
Остальные способы оптимизации, небольшие, изложу по мере раскрытия работы программы. Часть используемых техник взято из статьи Недокументированные возможности ПМК МК-61.
Как и в исходной программе результат хранится в регистрах R0…Re. Начало (R0) продублировано на экране (RX). В сумме 94 цифры, потому что в Re уже восемь цифр.
Для повтора вычислений нажмите В/0 С/П.
| # | | 00 | 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 |
|---|---|---|---|---|---|---|---|---|---|---|
| 00 | | В/О | 6 | 5 | F1/x | x→П0 | Cx | Кx→П0 | П→x0 | Fx=0 | 05 |
| 10 | | FВx | 7 | F10x | Fπ | КП→xe | ПП | E1 | 1 | + | x→Пe |
| 20 | | КП→xe | − | Fx=0 | 14 | П→xe | ВП | /-/ | 2 | К[x] | <-> |
| 30 | | F↻ | + | <-> | В↑ | FВx | <-> | ÷ | FВx | F↻ | К[x] |
| 40 | | Кx→Пe | × | − | 6 | F10x | × | П→xe | ВП | 2 | x→Пe |
| 50 | | + | КЗН | − | Fx=0 | 11 | П→xe | 2 | 5 | + | x→Пe |
| 60 | | 1 | Fex | x→П0 | С/П | ||||||
На входе в стеке предполагается:
| Регистр стека | Значение |
|---|---|
| RX | Очередные 6 разрядов для деления, обозначим как K. |
| RY | Игнорируем. Мусор. Обозначим как M. |
| RZ | Накопленный остаток, может быть 8 разрядный. Обозначим как S. |
| RT | Число N, на которое делим. Вначале 65. |
На выходе получается:
| Регистр стека | Значение |
|---|---|
| RX | Содержимое регистра Re. |
| RY | Новый S, как остаток от деления выше, уже умноженный на 106. |
| RZ | N. |
| RT | N. |
При этом, новое K = [(K+S)/N], результат целочисленного деления, уже будет сохранён в регистре, на которые указывает Re.
Кстати, небольшое уточнения. Делитель состоит из 2 знаков, остаток может тоже. При умножении на 106 уже 8 знаков. Это предел для ПМК. Поэтому группы по 6 знаков, а не по 7. В первой группе начинаем с 107, поэтому её это ограничение не касается. В первой группе может быть, а в конце и будет, 7 знаков.
Итак, стек выглядит как K, M, S, N. Теперь по порядку, что выполняем в подпрограмме:
| Действия | Стек | Адреса |
|---|---|---|
| Удаляем M, путём откидывания в хвост стека | K, S, N, M | 29..30 |
| K+S меняем местами с N | N, K+S, M, M | 31..32 |
| Дублируем первых два значения | K+S, N, N, K+S | 33..34 |
| Производим деление | (K+S)/N, N, K+S, K+S | 35..36 |
| Сохраняем N в конец стека | (K+S)/N, N, K+S, N | 37..38 |
| Получаем результат [(K+S)/N], и сохраняем | [(K+S)/N], N, K+S, N | 39..40 |
| Находим остаток S (новое) | (K+S) - ([(K+S)/N] × N), N, N, N | 41..42 |
| Умножаем остаток на 106 | S (новое) × 106, N, N, N | 43..45 |
| Извлекаем Re и пока не понятная зачем-то команда ВП | Re, S (новое) × 106, N, N | 46..47 |
И после адреса 47 код не похож на конец подпрограммы. Но…тут используем один из трюков, связанный с нелинейной адресацией при больших адресах. (см. Программное адресное пространство). По этой причине вызываем подпрограмму (см. адреса 15..16), как бы с адреса E1, который соответствует адресу 29. Но когда выполнение дойдёт до F9 (47), то автоматически произойдёт переход на адрес 00, где и располагается команда В/0. Рассмотрев подпрограмму, изучим всё от начала.
В/0 служит не только, как конец подпрограммы из абзаца выше, но и чтобы спокойно ввести N (=65), не опасаясь, что продолжится ввод пользователя.
Теперь нужно очистить 15 регистров. Это делают команды по адресам 05..09, но в R0 будет не 15, а 1/65 = 1.5384615^-02, которое и содержит 15 (в конце). Для косвенной адресации это не важно (см. Число положительное, но меньше единицы). Цикл очистки прервется, когда занулится сам R0.
Теперь командой по адресу 10 возвращаем 65, которое всё это время сохранялось в X1 и переходим к…
С адреса 11. Как упоминалось, Re используется как счетчик внутреннего цикла по всем регистрам. В самом начале это ноль. В RX хранится N.
Для подготовки вызова подпрограммы стек заполняется (адреса 11..14) значениями R0 (косвенно через Re), π, 107, N.
Тут мусор
в виде числа π вводим специально, чтобы соответствовать контракту подпрограммы.
Кроме того, такое значение удобно для отладки. Можно указать любую цифру или команду В↑.
После вызова подпрограммы, помня что в RX=Re, сразу прибавляем единицу и сохраняем новое Re (16..19).
Чтобы определить, что внутренний цикл закончен, косвенно извлекаем по Re. Если это само Re и есть, то дошли до конца. Что и делает проверка по адресам 21..23. Если нет, то с адреса 14 внутренний цикл повторится, только вместо π и 107, будут вот эта неудачная разность и новый S.
Внутренний цикл завершён. Значение Re обрабатываем отдельно. У Re нужно отрезать две последние цифры, что делается по адресам 24..28. Здесь и далее сознательно используем опасную в программном режиме ВП, вместо F10x, чтобы не сдвигать стек и не потерять значение N. К тому же и короче. После П→xe команда ВП безопасна (см. X2-влияющие команды).
А затем (адрес 29) подходим к подпрограмме, но напрямую. Она сделает, что должна, обновит Re. Но пройдём её напрямую, без возвратов. Тогда становятся понятными команды по адресам 46..49 – умножение Re на 100 и сохранение, для подготовки к следующему внутреннему циклу.
Теперь нужно уменьшить N, и проверить, что он не ноль. Сейчас в стеке: Re×100, S, N, N. Т. к. S может быть и нулём, то используем команду +, чтобы подтянуть стек и получить точно не ноль. Последний превращаем через КЗН в единицу и вычитаем из N, с проверкой по адресам 53, 54.
Тут довольно просто. К числу в Re, который уже умножен на 100, прибавляется 25 (эээ, ну просто знаю эти цифры числа e). Вместо длинного получения начала числа e из R0 (17182818), немного жульничаю и сразу вычисляю 2.7182818.
Всё. Почти. Стоит обратить внимание, что при первом проходе, когда Re ещё ноль,
после КП→xe
в Re станет не ноль, а 00000000. . Поэтому, после команды
ВП
по адресу 47, получится не ноль, а 10000000. .
Но… 91 это гарантированная точность алгоритма. В действительности она выше и есть запас прочности
.
Поэтому небольшое отклонение в самом начале, при вычислении минимальной прибавки и
для самых не значащих разрядов допустимо. Это не влияет на результат.